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ABSTRACT 

Context. We analyse the presence of nonradial oscillations in Cepheids, a problem which has not been theoretically 
revised since the work of Dziembowsky (1977) and Osaki (1977). Our analysis is motivated by a work of Moskalik et 
al. (2004) which reports the detection of low amplitude periodicities in a few Cepheids of the large Magellanic cloud. 
These newly discovered periodicities were interpreted as nonradial modes. 

Aims. Based on linear nonadiabatic stability analysis, our goal is to reanalyse the presence and stability of nonradial 
modes, taking into account improvement in the main input physics required for the modelling of Cepheids. 
Methods. We compare the results obtained from two different numerical methods used to solve the set of differential 
equations: a matrix method and the Ricatti method. 

Results. We show the limitation of the matrix method to find low order p-modes (I < 6), because of their dual character 
in evolved stars such as Cepheids. For higher order p-modes, we find an excellent agreement between the two methods. 
Conclusions. No nonradial instability is found below I — 5, whereas many unstable nonradial modes exist for higher 
orders. We also find that nonradial modes remain unstable, even at hotter effective temperatures than the blue edge of 
the Cepheid instability strip, where no radial pulsations are expected. 
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1. Introduction 

Cepheids are well known radial pulsators, oscillating in 
\ the fundamental mode or the low overtones. Despite a 
very abundant literature analysing theoretically and ob- 
servationally their radial pulsation properties, only very 
few studies have been devoted to nonradial oscillations in 
■ these stars. Dziembowski (1971) was the first to study non- 
radial oscillations in Cepheids, based on a quasiadiabatic 
linear stability analysis, and found that low-order modes 
(I < 2) are stable. Following this study, Osaki (1977) and 
Dziembowski (1977) showed that higher degree nonradial 
modes (I > 4) can become unstable. These analysis high- 
lighted the difficulty to calculate low-order p-modes in such 
evolved stars. These modes show a dual character as they 
can also propagate in the gravity-wave region, behaving 
as high-order gravity modes in the central regions and thus 
having several thousand nodes there. Since no observational 
evidence for the presence of nonradial modes in Cepheids 
was available at the time of these studies, no further the- 
oretical analysis of the stability of nonradial modes was 
performed. 

Later on, however, due to advances in spectroscopy 
and radial velocity measurement techniques, the presence 
of nonradial modes in classical Cepheids was more con- 
vincingly considered to explain anomalous behaviors such 
as amplitude modulations (Hatzes & Cochran 1995; Van 
Hoolst & Waelkens 1995; Koen 2001; Kovtyukh et al. 2003). 
These have been interpreted either as the beating of two 
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linearly unstable radial and nonradial modes with similar 
frequencies, or as a nonlinear interaction between a linearly 
unstable radial mode and a linearly stable, low degree, low- 
order nonradial mode (Van Hoolst & Waelkens 1995; Van 
Hoolst ct al. 1998). The latter hypothesis was theoretically 
investigated by Van Hoolst & Waelkens (1995), based on 
frequencies and linear growth rates of the nonradial modes 
calculated previously by Osaki (1977). 

More recently, Moskalik et al. (2004) found low ampli- 
tude secondary periodicities in a few first overtone Cepheids 
of the Large Magellanic Cloud, with frequencies close to the 
unstable, radial, first overtone. These newly discovered pe- 
riodicities were interpreted as nonradial modes. They pro- 
vide, to our knowledge, the first direct evidence for the 
presence of nonradial modes in classical Cepheids. But un- 
fortunately, these observations have not been confirmed yet. 
The observations of Moskalik et al. (2004) are based on pho- 
tometrical data, implying that if nonradial modes are in- 
deed detected, they must be low-degree modes. According 
to Osaki (1977), however, such modes should be stable. 
Motivated by the observations of Moskalik et al. (2004) and 
the fact that since the work of Osaki (1977), no updated 
stability analysis of nonradial modes has been performed, 
we reinvestigate this problem, taking into account the im- 
provement in the main input physics (opacities, equation of 
state) required for the modelling of Cepheids. 

We have used two different numerical methods for solv- 
ing the set of linear pulsation equations. The first one is 
based on a Henyey-type relaxation scheme, which is the 
most commonly used in the community; the second one 
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is based on the Riccati method (see Gautschy Sz Glatzel, 
1990, and references therein). The two methods are briefly 
described in §2 and §3 respectively. Because of the above- 
mentioned dual character of low-degree p-modes in highly 
evolved stars, the first method requires an extremely large 
number of grid-points for the Cepheid model because of the 
large number of nodes in the central region. Such method 
thus suffers from resolution problems and a lack of accu- 
racy. One of our purposes is to analyse and highlight such 
limitations and is done in §2. To overcome the problems, 
several ideas have been suggested, such as asymptotic meth- 
ods providing analytical solutions (Dziembowski 1977; Lee 
1985). Osaki (1977) found that the stellar envelope could 
be regarded as a unique pulsating unit, with an appropriate 
boundary condition at the bottom of the envelope which 
in principle, should take into account the enery leakage 
into the core. As a more attractive alternative, the Riccati 
method as used by Gautschy & Glatzel (1990), avoids the 
extra-complexity of asymptotic methods and can take into 
account the entire stellar structure, providing a correct de- 
scription of the energy leakage into the core without resort- 
ing to any approximation. 

In section §4 we present the outcome of our stability 
analysis and compare the results obtained with both nu- 
merical methods for low-degree p-modes. 

2. Linear stability analysis based on the 
Henyey-method 

2.1. Linear pulsation equations 

The pulsation calculations are performed with a nonradial 
code originally developed by Lee (1985) and previously used 
by one of the authors for an extensive study of radial pulsa- 
tions in classical Cepheids (Alibert et al. 1999). We briefly 
describe the method below. The system of linearised equa- 
tions describing small amplitude stellar pulsations can be 
found in Unno et al. (1989). The eigenfunctions X(r, 9, <j), t), 
solutions of this system of equations, are expressed in terms 
of spherical harmonics as X = x(r)Yj n (9, <f))e l ' Tt , where I in- 
dicates the degree of the eigenmode and m its azimuthal 
number, a = ay + ioy is the eigenfrequency of the mode, 
with P — 27r/<7 r the pulsation period and o$ character- 
ising the stability of the eigenmode. Positive values of Oi 
indicate stable modes. In the following, all frequencies are 
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, with R the radius and m 



the mass of the star. Convection is frozen in, assuming that 
the perturbation of the convective flux in the linearised en- 
ergy conservation equation is neglected. Though crude, this 
approximation can be justified as we will be mainly inter- 
ested in Cepheids close to the blue edge of the instabil- 
ity strip, where convection plays a minor role. As shown 
in Alibert et al. (1999), frozen-in convection provides a 
theoretical blue edge in good agreement with observed ra- 
dial fundamental and first overtone Cepheid pulsators. The 
boundary conditions are those imposed by the regularity of 
the eigenfunctions (see Unno et al. 1989). A normalisation 
condition is added : the radial component of the displace- 
ment is set to one at the surface of the star. We adopt, 

as in Unno et al. (1989), the following variables yi — — , 
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with £ r the radial and £/, the hor- 



izontal component of the displacement respectively, g the 
gravitational acceleration, <fi the gravitational potential, p' 
and 4>' the Eulerian perturbations of the pressure and grav- 
itational potential, L ra d and <5L ra d the radiative luminosity 
and its Lagrangian perturbation, SS the Lagrangian pertur- 
bation of the entropy, and C p the specific heat at constant 
pressure. The inner boundary conditions are: 
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and 



^2/3 -2/4 = 0, 
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(2) 



with ci = — -fej and m r is the mass inside the sphere of 
radius r. The outer boundary conditions at r=R are: 



(I + 1)2/3 + 2/4 = 0, (3) 
(2 - 4V ad U)2/i + 4V ad U(j/ 2 - y 3 ) + Ay 5 - y e = 0, (4) 
2A = L (5) 

with V a d the adiabatic gradient and V = — — . Finally, 

d In r 

the system of linearised equations is numerically solved with 
a Henyey-type relaxation method. It requires a guess for 
the eigenfrequency which is derived from the solution of 
the adiabatic problem. 

2.2. Modal classification 

The standard modal classification of nonradial modes, 
based on the determination of the number of nodes N g in 
the gravity-wave and N p in the acoustic- wave zones respec- 
tively (see e.g Unno et al. 1989) usually fails for evolved 
stars, as underlined by Dziembowski (1971), because of the 
dual character of the modes. Modes with p-mode charac- 
ter in the envelope do rapidly oscillate, like g-modes, in 
the central region. This behavior is illustrated in Fig. Q] 
which displays the radial displacement of a p-mode of de- 
gree Z=10 found in a typical Cepheid model. In this case 
the number of nodes calculated with the present numerical 
method is inaccurate and meaningless, with N p =117 and 
A r g =143. Another method must be adopted to select and 
classify p-modes, which are the most interesting ones since 
they can propagate to the surface and are thus potentially 
detectable. 

Also, a huge number of solutions is found, among which 
some are unphysical and result from numerical noise. For 
the selection of physical p-modes, we have used a criterion 
based on the modal kinetic energy, E c , which reads for a 
layer between radii r and r + dr: 

£c = ^VW, 

with|<5r| 2 = |e r | 2 + Za + l)ia| 2 - 

If the kinetic energy of an eigenmode reaches its max- 
imum value in the acoustic-wave zone (p > N,L, with N 
and L being the Brunt- Vaisala and the Lamb frequencies 
respectively), it is selected as a correct solution, which is 
illustrated in Fig. [21 Each selected eigenmode, for a given I, 
1 dip' is then classified as fundamental mode F or overtone (1H, 
^ ' 2H, etc) with increasing value of ay. This selection method 
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Fig. 1. Radial displacement as a function of temperature 
of a p-mode with Z=10, frequency oy = 3.62 in units of Co, 
for a typical Cepheid model (to = 5M Q , log L/Lq = 3.1, 
T off = 5930 K). 



allows us to find p-modes with high-degree, I > 6, but fails 
for lower degrees due to the limitation of the numerical 
method. Although we find numerous solutions for I < 6, 
many of them have very close values of a r but very dif- 
ferent values of cr^. Moreover, they all have pathological 
eigenfunctions and a selection of the correct solution based 
on the shape of the eigenfunctions or the kinetic energy is 
not possible. 

Since the problems above mentioned stem from the 
rapid oscillations of the eigenfunctions close to the cen- 
ter of the star, one could think of using the quite common 
"envelope model" approach, namely solving the dispersion 
equation using envelope models rather than complete stel- 
lar models. This approach, however, is not useful in the 
present context since one needs to transfer the unambiguous 
inner boundary conditions taken at r = to the inner edge 
of the envelope. Since there is no zone where the evanescent 
oscillation would be rapidly switched off, it is quite difficult 
to choose the correct inner boundary conditions with enve- 
lope models. In order to get rid of this problem, the inner 
boundary conditions should be taken deep enough, if not, 
they will affect the solution. We found that using enve- 
lope models requires the same typical resolution as the one 
required with complete stellar models. There is thus no ad- 
vantage of using an envelope model rather than a complete 
one since the size of the matrix which needs to be inverted 



Fig. 2. Kinetic energy E c (solid line, arbitrary units) for 
the same eigenmode with /=10 as in Fig. [T] E c reaches its 
maximum value in the p-zone, where the eigenfrequency 
(dotted line) is greater than the Brunt- Vaisala frequency 
N (dashed line) and greater than the Lamb frequency L 
(dash-dotted line). All frequencies are in s _1 (right-hand 
y-axis). This eigenmode fullfills the selection criterion of a 
p-mode (see text §2.2). 

is approximately the same. If the spatial resolution of en- 
velope models is too low, this implies adopting arbitrary 
inner boundary conditions, leading to strongly uncertain 
results for the values and the sign of <ii. The analysis of 
non-radial pulsations in Cepheids based on envelope mod- 
els should thus be used with caution. Finally, note that all 
results presented in this paper with the Henyey-method are 
obtained with complete stellar models and a typical num- 
ber of grid-points of 5000, in agreement with the required 
resolution estimate based on a WKB approach. 

3. The Riccati method 

Because of the difficulties encountered with the above- 
mentioned numerical method, an independent check of the 
results both for high and low degree p-modes is mandatory. 
Thereby it will be of particular interest to check whether 
the solutions selected with the previously described ma- 
trix method are correct and no physical solution has been 
missed. 

A method which does not suffer from resolution prob- 
lems and provides eigenvalues and eigenfunctions with pre- 
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scribed accuracy irrespective of the order of an eigenfunc- 
tion is the Riccati method. Its application to stellar sta- 
bility problems is described in Gautschy & Glatzel (1990) 
for radial perturbations and in Glatzel & Gautschy (1992) 
for nonradial perturbations. We shall only briefly summa- 
rize the basic essentials of the approach here and refer the 
reader to these publications for details. 

The rapidly oscillating character in the central regions 
of the eigensolutions considered here implies severe reso- 
lution problems due to the limited size of matrices, if the 
eigenvalue problem is solved on the basis of a matrix eigen- 
value problem. This basic difficulty, which always occurs 
for high order eigenfunctions can in principle be overcome 
by shooting methods, which do not rely on the inversion 
of large matrices. Moreover, using a shooting method, the 
local stepsize can be adjusted to match prescribed accuracy 
requirements. However, for higher order boundary value 
problems (the nonadiabatic nonradial stability analysis im- 
plies the solution of a sixth order complex differential equa- 
tion) shooting methods tend to suffer from the parasitic 
growth problem and are numerically unstable. Moreover, 
with increasing order the number of parameters which need 
to be iterated in the shooting process becomes prohibitively 
large. Therefore shooting methods are usually not consid- 
ered for the solution of high order boundary value problems. 

One way to overcome the parasitic growth problem 
and to reduce the parameters to be iterated to the com- 
plex eigenvalue while still preserving the basic advantages 
of a shooting method consists of transforming the linear 
boundary eigenvalue problem to a nonlinear one imply- 
ing a Riccati type equation. For the shooting process we 
thus obtain unambiguous initial conditions at both bound- 
aries and the only parameter which needs to be iterated 
is the eigenfrequency. In addition, integration of the non- 
linear Riccati equation as an initial value problem is nu- 
merically stable. As a result, the integration of the Riccati 
equation provides a complex (determinant) function of the 
complex frequency. Zeros of this determinant correspond to 
the eigenvalues sought. 

The latter property of the Riccati method offers an addi- 
tional advantage: eigensolutions can be determined without 
resorting to any initial guess. Using matrix methods, these 
are usually obtained on the basis of approximations (e.g., 
the solution of the adiabatic problem) which may lead to a 
missing of unexpected solutions. Using the Riccati proce- 
dure, always the full set of equations without any approx- 
imation is integrated to provide the determinant function. 
By tabulating the determinant as a function of the complex 
frequency, a coarse determination of the eigensolutions can 
be done which subsequently may be improved by iteration. 
Due to the existence of the determinant function there is 
no problem of spurious eigenvalues when using the Riccati 
method. 



4. Results 

4.1. Comparison of the results between the two methods 

We have compared results obtained with the two methods 
described respectively in §2 and §3 for a specific Cepheid 
model. Since our work is motivated by the observations of 
Moskalik et al. (2004), we selected a model which could 
describe the properties of one of the observed LMC tar- 
gets where secondary periodicities were detected. The ob- 



Table 1. Linear stability analysis results for a Cepheid 
model with mass m = 5M Q , metallicity (in mass fraction) 
Z = 0.01, log L/Lq = 3.1, T cff = 5930 K. The real part 
ay and imaginary part Oi of the eigenfrequencies are nor- 
malised to do ~ 1.42 10 5 cgs. Results obtained with the 
Riccati (Ric.) and the Henyey (Heny) methods are given. 
Positive values of Ui indicate stable modes. 



1 mode 


Or 


Oi 


Or 


Oi 


Period 




(Ric.) 


(Ric.) 


(Heny.) 


(Heny.) 


(days) 


F 


2.99 


-3.8 10" 4 


3.02 


-4.4 10~ 4 


3.43 


1H 


4.23 


-6.3 10" 3 


4.25 


-7.5 10~ 3 


2.43 


2H 


5.40 


-2.7 10~ 3 


5.40 


-4.2 10~ 3 


1.92 


6 F 


3.13 


2.0 10~ 2 


3.20 


1.8 10~ 2 


3.24 


1H 


4.57 


-3.9 10" 3 


4.58 


-5.0 10" 4 


2.26 


2H 


5.98 


-1.1 10~ 3 


5.96 


7.5 10" 3 


1.74 


10 F 


3.60 


-2.5 10~ 3 


3.61 


-1.9 10~ 3 


2.87 


1H 


5.23 


-9.8 10~ 3 


5.22 


-5.2 10~ 3 


1.98 


2H 


6.88 


2.4 10" 2 


6.83 


3.8 10" 2 


1.52 





served Cepheid, LMC-SC2-208897, is a first overtone pul- 
sator with a period of Pi=2A2 days and absolute magni- 
tude My ~ —3.2. Adopting the same input physics and 
evolutionary code as in Alibert et al. (1999), we found that 
a model with mass m = 5Mq, metallicity (in mass frac- 
tion) Z = 0.01, log L/Lq = 3.1, T cff = 5930 K and close to 
the end of central He burning (central mass fraction of He 
Y c ~ 4. 10 ~ 2 ) provides a good fit for the observed magni- 
tude and period of LMC-SC2-208897. The calculated first 
overtone period is unstable with a period Pi =2.43 days and 
the absolute magnitude is My = -3, in excellent agreement 
with the observed values (see Alibert et al. 1999 for the 
determination of absolute magnitudes). 

Excellent agreement is found between the two meth- 
ods for the periods and growth rates of radial modes (see 
Table 1). For nonradial modes with degree I > 6, Fig. [3] 
compares the values of oy and <jj obtained with both meth- 
ods, up to 1—20. The values of ay agree within less than 
2%. Differences are found for the values of cr, for over- 
tones, whereas the agreement is excellent for the funda- 
mental mode (see Fig. [3]). Values of o r and cr, for a few 
nonradial modes are also given in Table 1. 

For the above-mentioned effective temperature and lu- 
minosity, a variation of the mixing length, between one 
pressure scale height Hp and 2 x Hp or of the metallicity 
between Z=0.0l and Z=0.02 have only a slight influence: 
OY changes by less than 2.5% and <Ji by less than 8 % (if 
\<Ji\ > 5.10 -3 ). Larger differences can appear, however, for 

very small values of Oi (\<Ji \ <, 5.10 -3 ), when using the ma- 
trix method described in §2. 

We have also compared the results obtained with both 
methods for different effective temperatures, adopting the 
same mass and luminosity as the above-mentioned Cepheid 
model. Fig. [5] shows the results for the fundamental mode, 
1H and 2H with degree 1=10. Here again, excellent agree- 
ment is found for the values of oy , but differences appear for 
the values of Oi for the overtones. In general, both methods 
agree on the sign of <7j, providing the same results concern- 
ing the stability properties of modes, except for very small 
values of <n (lovl < 5.10~ 3 ). 
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Fig. 3. Comparison, for our selected Cepheid model (m = 
5M , log L/L Q =3.1, = 5930 K) between the results 
obtained with the Henyey method (see §2, dotted lines) and 
the Riccati method (see §3, solid lines). The upper panel 
displays the real part of eigenfrequencies of p-modes as a 
function of the degree I and the lower panel the imaginary 
part, in units of ao- Each curve is labelled, on the right hand 
side, according to the modal classification. Note that for the 
fundamental mode, the curves corresponding to the values 
of Oi obtained with both methods are indistinguishable. 



Both for low and high degree modes we have tested 
the dependence on the outer boundary conditions of the 
results of the stability analysis. This is necessary, since 
the outer boundary of the stellar model (the photosphere) 
does not correspond to the physical boundary of the star. 
Therefore the thermal and mechanical boundary conditions 



Fig. 4. Variation of ay and en, in units of ao, as a func- 
tion of effective temperature for models with m = 5M , 
logL/L Q = 3.1 and Z=0.01. Different p-modes with degree 
1= 10 are displayed, as indicated on the figure. The solid 
line corresponds to the Riccati method and the dashed line 
to the Henyey method. 



are ambiguous there. In fact, both for radial and nonradial 
modes we found a sensitive dependence on boundary condi- 
tions even within the physically admissible ones. The results 
shown in the paper, however, are based on the most conser- 
vative condition with respect to stability. I.e., we have cho- 
sen the condition which provides the least unstable modes. 
It corresponds to the requirement of vanishing Lagrangian 
pressure and temperature perturbation. All other boundary 
conditions, in particular at the inner boundary, are unam- 
biguous. 
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In contrast to the thermal and mechanical outer bound- 
ary conditions the influence of metallicity on the stability of 
all the models investigated is rather weak. Relative differ- 
ences between the results for Z=0.0l and Z=0.02 amount 
to 10 per cent at maximum for Oi. 

4.2. Results for low-degree p-modes 

Modes of degree below I w 6 are extremely difficult to find 
with matrix methods. Indeed, only high overtones can be 
found (e.g. 5H for I = 2). Some results obtained in this 
regime are based on the Riccati method and shown in Fig. 
[5] . Both the frequencies of p- and g-modes decrease with 
the harmonic degree. However, while p-modes reach finite 
frequency values for I — * 0, the frequencies of g-modes van- 
ish in this limit. The stellar models considered here exhibit 
a propagation region for gravity waves close to the center 
of the star which for high values of the harmonic degree 
is well detached and shielded from the acoustic propaga- 
tion region by an efficient evanescent barrier. It allows for 
gravity modes with frequencies in the range of p-modes up 
to quite high order. Since the decrease of frequencies with 
I is stronger for the g-modes, this leads to multiple cross- 
ings and resonances between p- and g-modes. Due to the 
efficient evanescent barrier for high values of I the interac- 
tion of p- and g-modes is weak at resonances above I « 6 
and leaves them unaffected. For I < 6, however, the barrier 
becomes weaker and the resonances at the crossing of g- 
and p-modes imply significant interaction and unfold into 
avoided crossings. Bumps both in the real (less pronounced) 
and imaginary parts of the eigenfrequencies as a function 
of I found in Fig. [5] are parts of these avoided crossings. For 
illustration a part of the run of a g-mode with I is shown 
there as a dashed line and I is taken as a real parame- 
ter. The latter allows one to follow modes continuously (of 
course, only integer values are physically meaningful). As a 
consequence of the resonances, modes with I < 6, in partic- 
ular the so-called p-modes, should not be classified as either 
p- or g-modes. Physically they rather exhibit properties of 
both types of modes. 

5. Discussion and conclusion 

The stellar model whose parameters are thought to match 
the properties of the Cepheid observed by Moskalik et al. 
(2004) is found unstable both with respect to radial and 
nonradial perturbations. Following modes up to I = 500, 
we find the following unstable modes: the fundamental p- 
mode for I = and between 9 < Z < 368, the first overtone 
for I = and between 6 < I < 13 and the second overtone 
for I =0, 5 and 6. No nonradial instability was found below 
/ = 5: Specific studies based on the Ricatti method for I = 1 
and 2 have not revealed any instability in the p-mode range. 
Moreover, from extrapolating the curves in Fig. [5] we do not 
expect instability for I = 3 and 4. We can thus conclude that 
if low degree p-modes with I < 2 were necessary to explain 
the observations of Moskalik et al. (2004), they cannot be 
attributed to a linear instability of these modes. Rather 
other effects, such as, e.g., nonlinear mode coupling, have 
to be invoked for an explanation of their excitation. 

With respect to the dependence on effective tempera- 
ture of the instability of radial and nonradial modes we find 
the unstable radial modes (F, 1H, 2H) to stabilize above 
6200K, for our test case luminosity log L/L & — 3.1. The 
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Fig. 5. Variation of a r and at, in units of do, for the modes 
indicated as a function of the harmonic degree for low values 
of I and the Cepheid model having m = 5Mq, logL/L Q = 
3.1, T cff = 5930 K and Z=0.01. All results are based on the 
Riccati method. The bumps for low values of I are due to 
an interaction of p- and g-modes through avoided crossings. 
The dashed line indicates the run of a selected g-mode. 

same is true for nonradial modes with I < 10. However, 
fundamental p-modes with I > 10 may still be unstable at 
temperatures above 6200 K, where no radial mode is found 
to be unstable. For example, the fundamental p-mode for 
Z = 15 is unstable up to 6290 K, the corresponding mode 
for Z = 60 is unstable up to 6640 K. Although observation- 
ally difficult to detect, we thus conclude that it might be 
worthwhile to search for nonradial pulsations with / > 5 in 
Cepheids. Particularly promising appear to be observations 
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of objects with T e g > 6200 K, where no radial pulsations 
are expected and - according to theory - only nonradial 
pulsations should prevail. Observations of nonradial pul- 
sations in Cepheids could then provide an interesting and 
unique tool to sound their inner structure. 
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